Randomness-Induced Redistribution of Vibrational Frequencies in Amorphous Solids 
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Much of the discussion in the literature of the low frequency part of the density of states of 
amorphous solids was dominated for years by comparing measured or simulated density of states 
to the classical Debye model. Since this model is hardly appropriate for the materials at hand, 
this created some amount of confusion regarding the existence and universality of the so- called 
"Boson Peak" which results from such comparisons. We propose that one should pay attention to 
the different roles played by different aspects of disorder, the first being disorder in the interaction 
strengths, the second positional disorder, and the third coordination disorder. These have different 
effects on the low-frequency part of the density of states. We examine the density of states of a 
number of tractable models in one and two dimensions, and reach a clearer picture of the softening 
and redistribution of frequencies in such materials. We discuss the effects of disorder on the elastic 
moduli and the relation of the latter to frequency softening, reaching the final conclusion that the 
Boson peak is not universal at all. 



I. INTRODUCTION 

The study of the density of states of solid materials 
started with attempts to understand the temperature de- 
pendence of the specific heat at low temperatures, say 
Cv = (dU/dT)v where U is the energy and T the tem- 
perature of the system. This called for a microscopic 
theory for solids, and the first one was developed by Ein- 
stein, assuming that in d dimensions each atom is rep- 
resented as a d-dimensional harmonic oscillator [l[ (in 
the original paper the case d = 3 was considered). In 
this article Planck's quantization assumption, which was 
originally applied to radiation, was extended to solid vi- 
brations [2| . In the case of dN linear oscillators each with 
its own frequency W{, Einstein's result can be expressed 
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Here ks and K are Boltzmann's and Planck's constants 
respectively, and the density of states g(uj) is defined by 
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where b(x) is the delta function. 

Both the theoretical calculation and experimental mea- 
surement of g(co) attracted enormous attention over the 
last century. We are interested in amorphous system 
like glasses, gels, foams etc., in which randomness ap- 
pears to influence the low-frequency part of the density 
of states g(u>). In particular the relation between the low- 
frequency behavior and the low-temperature thermody- 
namics of such systems is of great interest. Studies of the 
low-frequency part of the density of states are dominated 
by dividing g(u>) by the prediction of the Debye model, 



focusing on the deviation between the two, and in par- 
ticular on the so-called "Boson peak" which emerges in 
many cases. There exist numerous claims about the Bo- 
son peak, its universality Q and its relation to softening 
or hardening of the materials under changes of material 
parameters [J]. In this paper we first explain the classi- 
cal approaches to the issue, including the Debye model 
and beyond, and then we examine the issue of universal- 
ity of the Boson peak in one and two dimensions using 
tractable models that can be computed to desired accu- 
racy. We conclude that there is nothing universal about 
the Boson peak, and that different types of disorder re- 
sult in very different redistributions of the low-frequency 
modes over the spectral domain. There is in general no 
correlation between the size or the position of the Boson 
peak and the increase or decrease of elastic moduli. 

The structure of this paper is as follows: In Section ITTl 
we review Debye's theory and the historical origins of the 
Boson peak. In Section IIIII we remind the reader what 
is entailed in computing the density of states by solving 
the appropriate eigenvalue equation. We then remind the 
reader that even for a perfect cubic crystal the Debye 
model is not exact, with corrections at frequencies which 
become lower as the material gets softer. To understand 
the effect of disorder in the inter-particle forces we review 
some known results and present some new results for one- 
dimensional chains in Scction|lVj In Section[V]we discuss 
tractable models of disorder in two dimensions, aiming 
to better model the typical disorder exhibited by glass- 
forming systems. We first examine the effects of disorder 
in the spring-constants or, equivalently, in the positions 
of the particles. Second, we consider disorder in the co- 
ordination numbers (the number of nearest neighbors), 
demonstrating that this can lead to major corrections to 
the Debye form and to very large Boson peaks. This is 
in general agreement with the idea that the scenario of 
glass-formation can be encoded by the changing coordi- 
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nation numbers as a function of temperature (so-called 
upscaling [1 i, 0, B i, El EH)- In Section ED the elas- 
tic moduli and the Debye frequencies which define the 
'Debye point' are discussed. In Section IVIII we offer a 
summary of the paper and some concluding remarks. 



II. DEBYE'S MODEL AND THE BOSON PEAK 

A. Debye's Model 

For the sake of simplicity in comparing with experi- 
mental results, Einstein employed a minimal model in 
which all the oscillators have the same frequency oje [H- 
Then g(uj) = 5{lu — oje) and the specific heat is given by: 



of longitudinal and transverse sound waves: 




(3) 



where 9e = Swb/^b is the so-called Einstein tempera- 
ture. This minimal model agreed qualitatively with ex- 
perimental observations of the specific heat; neverthe- 
less careful measurements show that the low tempera- 
ture behavior of ^ for a three-dimensional solid, i.e. 
CV ~ 3Nk B (8 E /T) 2 exp(-8 E /T) falls off faster than ex- 
perimental values. To explain the observed data one 
needs to account for more adequate vibrational spectra of 
actual solids. Debye was the first to link the oscillator fre- 
quencies in {1} with the collective vibrations of the solid. 
He treated a solid as an elastic isotropic continuum and 
estimated the density of vibrational states for the case 
of a spherical body [!§]■ Later it was shown that the 
frequency distribution is independent of the shape if the 
size of a body is large enough and the surface contribu- 
tion can be neglected [l3| . The modern simple derivation 
of the Debye distribution can be found, e.g., in jljj . The 
main statement of the continuum approximation is the 
linear dispersion relation between the frequency and the 
absolute value of the wave vector fc, to = u s k where u s is 
the sound velocity. In an isotropic body there exist one 
longitudinal and d— 1 transverse sound waves, therefore, 
for spatial dimension d > 1 one needs two dispersion re- 
lations. The Debye density of states is given by: 
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where lod is the Debye frequency which defines the cut-off 
frequency in the spectrum: 
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Here VL d = (7r d / 2 )/r(l + d/2) is the coefficient in the vol- 
ume definition of d-dimensional hypersphere of radius r, 
Vd = H^r , is the gamma function, p is the particle 
number density, ui and u t are the speeds of propagation 
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where K and p are the bulk and shear moduli respec- 
tively and to is the molecular mass. 

Debye's model has the advantage of a simple analyt- 
ical form depending on the elastic properties of a solid 
only. Due to the long- wavelength approximation it is in- 
sensitive to the microscopic structure. Nevertheless ex- 
perimental measurements indicated from the start that 
Debye's model is far from being the end of the story. 



B. The Boson Peak 

The vibrational properties of solids can be investigated 
experimentally by studying the inelastic interactions of 
external radiation with the solid vibrations. For inelastic 
scattering of photons one observes the Raman effect, dis- 
covered by Raman in liquids 15] and by Landsberg and 
Mandelstam in crystals [l6j]. As a result of this effect 
the frequency of the incident photon is either red shifted 
(Stokes scattering, with high amplitude) or blue shifted 
(anti-Stokes scattering, with low amplitude). 

In crystals, due to the periodic structure, selection 
rules give rise to a discrete set of lines. In amorphous 
materials these spectral lines broaden, giving rise to a 
continuous spectrum. The Raman line shape was related 
to the density of states of amorphous materials in [l7| 
under some assumptions in the harmonic approximation. 
The result can be rewritten for Stokes scattering in the 
following form: 



L e -x p 
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(7) 



where I exp (uj) is the observed Raman intensity at the 
frequency shift equal to u) and 



n(ui) 
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(8) 



is the Bosc distribution function. The function C(oj) is 
an empirical function called "the average light vibration 
coupling constant" . Thus the right-hand side of Eq. ([7]) 
is independent of temperature, meaning that the tem- 
perature dependence of the Raman intensity should be 
compensated by the temperature dependence of the Bose 
distribution function. This conclusion was confirmed 
by experiments (see e.g. [3). At low frequencies the 
Raman spectrum has a bump whose amplitude changes 
with temperature. Once scaled by the Bose function the 
data at different temperatures collapse to a temperature- 
independent peak, which is therefore usually referred to 
as the 'Boson peak'. 

Analysis of Raman spectra for different amorphous 
materials indicates the existence of a Boson peak (l9j . 
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Therefore, it was suggested that Raman spectra indicate 
some universal features of amorphous systems, indepen- 
dent of the details of molecular interactions. Unfortu- 
nately, experimental results determine only the product 
of the density of states and the light-vibrational cou- 
pling constant, cf. Eq. (J7J). Under the assumption 
that the low-frequency density of state is defined by 
the Debye model, for three-dimensional systems we have 
g(uj)/u! 2 = const. In this case, and only in this case, Eq. 
([7]) implies that the coefficient C(u>) must exhibit 
the Boson peak. On the other hand if the Debye model 
does not apply to the particular material at hand, this 
conclusion cannot be reached. 

Additional light was shed on this problem using inelas- 
tic scattering of cold neutrons [13j . Such experiments in- 
dicate that in amorphous solids at small frequencies the 
quantity g(u))/u> 2 is not constant, showing an excess in 
the vibrational density of states [lo, [Il[ . The vibrational 
density of states defines the temperature dependence of 
the specific heat (JT]) and indeed the latter quantity also 
displays an excess at low temperature amorphous solids 
compared to the prediction of Debye's model. This again 
provides evidence for additional contributions to the vi- 
brational density of states [22j ■ 

At present it is clear (see, e.g. 23]) that the Raman 
coupling coefficient C(u>) is a rather complicated mono- 
tonic function of frequency. Therefore, the peak in ([7]) 
is defined by the non-monotonic behavior of the density 
of vibrational states at low frequencies (in the sense of 
Raman scattering). The term 'Boson peak' is transferred 
from the Raman intensity to the shape of the vibrational 
density of states at low frequencies. In other words, the 
Boson peak describes the deviation (excess) from the ex- 
pected constancy of g(ui)/uj d ^ 1 in the <i-dimensional De- 
bye model. 

Once we define the problem of the Boson peak as equiv- 
alent to finding the deviations from the Debye model, the 
Boson peak is no longer special to amorphous solids. De- 
bye's model takes into account only homogeneous elastic 
effect; after all, it is well known that the vibrational spec- 
tra of crystals have maxima at the van Hove singularity 
points and these maxima are independent of the tempera- 
ture. The spectral properties in these regions are defined 
only by the lattice structure and the dimensionality (see, 
e.g., the well-known exact solution Eq. (T53|) below). 

Disorder brings about additional deviations from the 
Debye model and different kinds of disorder have differ- 
ent effects on the density of states. We will show below 
that in one-dimensional systems disorder of the inter- 
particle interactions induces a frequency redistribution 
with smoothing of the van Hove singularity. The peak of 
the spectrum moves to low frequencies with a shift which 
depends on the distribution of interactions. The same re- 
sults were obtained by direct solution of Eq. (TT5]) below 
for three-dimensional cubic lattices with spring constants 
distributed in accordance with a Gaussian [24j or other 
distributions [25| (in contrast to one-dimensional chains, 
three-dimensional cubic lattices are stable even if some 



of the spring constants are zero or negative). These and 
other results using the coherent potential approximation 
(26| lead to a conclusion that the 'Boson peak' in dis- 
ordered systems is associated with the lowest van Hove 
singularity in the spectrum of the reference crystal [27| . 

It is important to stress that all these results can be 
taken only as a general indication for the appearance of 
the 'Boson peak' in amorphous solids in two or three di- 
mensions. In all these models only nearest neighbor har- 
monic interactions (spring constants) were considered. 
For cubic lattices in two and three dimensions such inter- 
action cannot give rise to a shear modulus, and only the 
bulk modulus is non-zero. Next-nearest-neighbor inter- 
action are necessary for having a non-zero shear modulus. 



III. BEYOND THE DEBYE MODEL: A FAIR 
WARNING 

A more general microscopic model of vibrations in a 
solid was proposed by Born and von Karman [28| . In the 
frame of this model it is assumed that all the atoms in 
a crystal interact with spring-like forces and that they 
vibrate near fixed equilibrium positions. This harmonic 
approximation can be used both for crystals and amor- 
phous solids; however, an analytical solution can be ob- 
tained only for very simple cases of regular crystal struc- 
tures. 

The total potential energy of a particle configuration 
R = {f\, r*2, ■ ■ ■ , f/v} is expressed in a pairwise approxi- 
mation as a sum over pair potentials: 

Relative particle positions are given by vectors r*y = fj — 
fi , the distance between the ith and jth particles is = 

Small displacements of all particles fi — > r\ = fj + 
Sfi lead to a new configuration R' = {r'i,r'2, . . . 
with the total potential energy: 

Ur> = \Y,4>(\^ +6 ^ I). ( 10 ) 

where: 

6f%j = Sfj - 5n. (11) 
We use the Taylor expansion 

cp(\f+8f\) = <t){r) + {5r-S7)(j)(r) + ^(5r-V) 2 <j>(r) + . . . (12) 

to obtain 

(13) 
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where 



II , , ' ' ' •'./ 'J - ;; | 4>_ ( r ij ) J 



(14) 



is a symmetric tensor = Tji, riij — rij fry, and X is 
the identity tensor. 

Substitution of (jlip to (fl~3f yields the dependence the 
energy of a harmonic system on particle displacements: 



Ur< = E/fl - E*^* ' ^ + ^ E 5fj ■ Aj • ^ 



(15) 



Here the force applied to the ith particle is defined by: 

(16) 



ij J'Hj 



and the dynamical matrix is given by 



J2Tik , if i = j 
-%j , if i ^ j 



(17) 



The equations of motion for a harmonic solid follow 
from (fTol): 



d 2 5rf 



P ■ Srf, 



(18) 



where rtii is the mass of ith particle and 1 < a, < d. In 
equilibrium Tf — and these equations are simplified. 

Substitution of a particle displacement of the form 
5rf = ufexp(-iwt) reduces Eqs. (fTS|) to the eigenvalue 
problem: 



L0 2 uf 



a/3 /3 
- ■ • 7/ . 



(19) 



J,/3 



This equation can be solved directly by diagonalizing 
t)^ for a system of N particles when N is not too large. 
Binning the resulting eigenvalues leads to a histogram 
that approximates the density of states. The simplest ex- 
ample is the d-dimcnsional cubic lattice with unstressed 
distance a between adjacent lattice points at zero pres- 
sure. In the approximation of nearest-neighbor interac- 
tions with spring constants 4>"{a) = 7, the matrix (Ti~4")) is 
given by: 



/TftQ 



7 , if I l a - m 
, otherwise, 



1=1 



(20) 



where a particle position is defined by the d-dimensional 
vector la with components {l a a} where l a are integer 
numbers. For this case the density of states can be found 
analytically [13j in the form of an inverse Laplace trans- 
form, 



e w s F(s)ds, 



(21) 



where the image function F(s) is: 

F(s) = i e -^% d (27 S ). 



(22) 



Here 7 = 7/ to and /g is the modified Bessel function of 
order zero [291 ]. The inverse Laplace transform in closed 
analytical form is defined for one and two dimensional 
systems, for example if d = 1 the density of states is 
given by: 



, < U) < UJ n 
, U) > LO m ax 



(23) 



where Lo max — 2^/7. This function diverges at u> = 0J m ax- 
This is a general property of the density of the vibrational 
states; for periodic structures there are integrable singu- 
larities (called van Hove singularities) of g(uS) (d = 1,2) 
or its derivatives (d — 3). The positions and types of van 
Hove singularities depend on the spatial dimension and 
the topological properties of the crystal [HI, [3(| • The low 
frequency behavior of this density of states was computed 
in [13| with the final result 



d 2 r(d/2) (47r 7 ) d / 2 



(24) 



The comparison of (|24[) with the Debye result ^ shows 
that the Debye model gives the first term in a more gen- 
eral expansion. The softer the system is, the larger is the 
correction. Substitution of Eq. (|24|) to Eq. (TTJ) allows 
to estimate corrections to Debye's specific heat (see, e.g., 
[3l|). We note that even for a perfect cubic crystal the 
Debye model is not exact, and there can be significant de- 
viations. Clearly when the crystal is not perfect or when 
disorder sets in the changes from the Debye limit can be- 
come much larger. Thus a blind comparison of any given 
density of states to the Debye limit may be unwarranted 
and can lead to spurious conclusions. We will come back 
to this issue when we discuss the Boson peak below. 



IV. ONE-DIMENSIONAL DISORDERED 
CHAINS 

Consider a one-dimensional harmonic chain with lat- 
tice spacing a and with random masses and spring con- 
stants as the simplest model of a disordered solid. In the 
nearest-neighbor approximation the interaction potential 



,■) is defined by: 



) = 7:7M+i( r M+i - a Y 



(25) 



where 7^+1 are random spring constants taken from a 
prescribed distribution ^(7). In this case the definition 
(fT4|) reads: 



_ J 7y , if I i-j |= 1 
, otherwise 



(26) 
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and Eqs. (fT9|) are written as: 

m l w 2 w i = —ji-i.iUi-i + 

(li-l,i + ~ 7i,i+l u »+l- (27) 

Unfortunately, it is impossible to derive a dispersion re- 
lations from these equations and the analytical solution 
discussed above becomes meaningless. Nevertheless, the 
response of the system to an applied static force AP can 
be inferred from the equilibrium conditions which follow 
from (1251): 



7*,i+i ( r i,i+i - a) = AP. (28) 
Summing Eqs. (|28p yields the elongation of the chain: 



JV-l 

AL = 2J r M+i 



w-i 



(iV- l)a 



7i,i+l 



(29) 



It is suitable to introduce a quantity = (-), then the 
bulk modulus defined by the condition (|29p is given by: 

(30) 



K = 



lav 
P 



We reiterate that i av is the harmonic average of 7. Sub- 
stitution of PD|) to © and to (O yields the following 
Debye frequency: 




(31) 



If lav > the Debye frequency has a finite value. Since 
the Debye model takes into account only the elastic prop- 
erties of the material, it should be exact in the limit 
u; — s- independently of the detailed structure of the ma- 
terial. In this limit every material is an elastic medium. 
Thus we expect lim^^o ff( w ) = rf/wo in agreement with 
the general law ([4J. We refer to this limit as the "Debye 
point" . 

In the case j av = the low frequency behavior of the 
density of states depends on the properties of the proba- 
bility distribution function p(i) for spring constants [32| . 
If p(i) 7 —,Q — > const in contrast to the Debye model 
the density of states exhibits the singular behavior 
g(u>) u —>o ~ V '—Incj. The density of states in the whole 
frequency region was obtained by Dyson in [33[ analyti- 
cally for a particular distribution of the ratio of the spring 
constants to the masses. Dyson introduced a set of new 
constants {in} defined by: 



l2n-l 



ln,n+l 



l2r, 



ln,n+l 

m n+ i 



(32) 



and derived an analytic solution for the distribution 



Pn(l) 



r(n) 



-n-l £ -n^ 



(33) 
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FIG. 1: Color Online: Density of vibrational states for a one- 
dimensional crystal (continuous line) and for exponentially 
distributed random interaction strength (dashed line). 



The frequencies are measured in units of y/(i/m), and 
(7) = 1. For asymptotically large n, p n {l) — ► 5(i — 1) 
which corresponds to the crystal state and the solution 
coincides with flSS]) [33j]. 

The density of states for the a special case of the dis- 
tribution (|33[) with n = 1 (exponential distribution) is 
shown in Fig. [T] It is known that disorder leads to 
smoothing out any van Hove singularity [26[ . The Dyson 
solution shows that due to the smoothing out of the peak, 
states penetrate into the high-frequency region which 
is forbidden for the periodic structure. In the small- 
frequency regime this function diverges logarithmically 
in accordance with the general result of [32j. The fre- 
quencies are redistributed so that the zero-frequency sin- 
gularity is followed by a dip. Such behavior is completely 
different from that of the Debye model. 

Unfortunately, it is impossible to study the crossover 
from Debye to non-Debye behavior analytically. Nev- 
ertheless the density of states of one-dimensional disor- 
dered systems can be estimated with the help of the ef- 
ficient numerical method proposed in I34H (extension for 
higher dimensions is discussed in [35|,[36|]). This method 
allows to calculate the number of frequencies less than u) 
using properties of a Sturm sequence [37] . In the follow- 
ing we present calculations pertaining to chains of 10 7 
particles of identical mass m = 1; in order to compare 
different systems we enforced in all cases cases (7) = 1. 
The results are summarized as follows: 



A. Uniform distribution. 

The simplest distribution function (used also in (34j for 
chains of 10 3 particles) is the uniform distribution: 



Pu(l) 



A < 7 < 1 
, otherwise 



(34) 
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FIG. 2: Color online: Density of vibrational states for a one- 
dimensional chain with interactions distributed by the uni- 
form distribution (|34|) . Different symbols pertain to different 
values of the parameter A, see inset. 



FIG. 3: Color online: Density of vibrational states for a 
one- dimensional chain with interactions distributed by the 
Weibull distribution (136[l , Different symbols pertain to dif- 
ferent values of the parameter a, see inset. 



For this distribution: 



2A 



h 



,1+A- 

'1-A 



(35) 



If A — > p u (l) — > 5{l — 1) and the system is reduced 
to the homogeneous chain. If A — > 1 the spring constant 
lav — ► and one can expect the divergence of g(u)) at 
vanishing frequencies. 

The density of vibrational states for the uniform distri- 
bution is shown in Fig. [21 Upon increasing the parameter 
A the van Hove singularity at u> = 2 in reduced units is 
smoothed out and then splits into two peaks moving in 
opposite directions. When A approaches unity the num- 
ber of low frequency modes increases and a minimum at 
intermediate frequencies is formed. 



B. Weibull distribution. 



The Weibull distribution is defined by: 

a-l 




The mean is given by: 

( 7 ) = AT(1 + 1/a) 



(36) 



(37) 



and in order to obtain (7) = 1 the parameter A was set 
to: 



A = 



1 



r(l + l/a) 
The average spring constant is given by: 

asinfa/a) 

lav = • 



(38) 



(39) 



When the parameter a — > 00, Pw(l) — * $(l ~ 1) an d 
the chain becomes uniform. For a = 1 the Weibull dis- 
tribution degenerates to the exponential distribution and 
the results obtained in [33| are expected. The density of 
the vibrational states for the Weibull distribution with 
a > 1 is shown in Fig. [3] In these cases one peak ad- 
vances toward low frequencies, and in the vicinity of zero 
frequency another peak is developed. 



C. Inverse distribution. 

The exponential distribution of the logarithm of the 
spring constant was used in [251 ] for the investigation of 
the vibrations of a three dimensional disordered cubic 
lattice. This distribution is given by: 



Ps{l) 



j\- , if a < 7 < Aa 
, otherwise 



(40) 



The mean is a(A — l)/ln\, therefore, the parameter a 
is defined by: 



InX 
a = . 

A-l 

The average spring constant is given by: 

2 



lav = A 



ln\ 
A-l 



(41) 



(42) 



A = 1 corresponds to the ordered chain. The densities 
of the vibrational states for the distribution (|4"01) with 
different values of A are shown in Fig. [5] For large A 
a dip at low frequencies develops and is followed by a 
peak, both moving to smaller frequencies with increasing 
the parameter A. 

Finally, recall that the Debye behavior of a disordered 
chain is defined by its elastic properties (|31[) . i.e., by i av . 
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V. TWO-DIMENSIONAL SYSTEMS 

In an amorphous solid the distance between two par- 
ticles, their relative orientation, and the number of near- 
est and next-nearest neighbors of each particle all possess 
some randomness and thus all the terms in 7y are ran- 
dom and so are the elements of the dynamical matrix X>y 
in Eq. (fT9| . In order to study the effect of different types 
of disorder on the spectrum of the dynamical matrix we 
examine several models of elastic networks. 



A. Elastic triangular anti-ferromagnet 



FIG. 4: Color online: Density of vibrational states for a one- 
dimensional chain with interactions distributed by the inverse 
distribution (|40|l . Different symbols pertain to different values 
of the parameter A, see inset. 
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FIG. 5: Color online: Comparison of the vibrational density 
of states for a common value "/ a v = 0.6838. The symbols are 
explained in the inset. 



Therefore it is useful to compare vibrational properties of 
different chains with the same bulk modulus, cf. Fig. [5] 
Note that this figure is interesting from the point of view 
of comparing with Debye's model. The Debye point at 
lu = is the same for all three models since we chose "f av 
to be the same; hence the Debye frequency (j31"j) is iden- 
tical for these three models. We could therefore expect 
identical Debye predictions for these three models. In 
contrast, the actual density of states presents widely dif- 
ferent frequency dependence for the three models. This 
means that the re-distribution of frequencies depends on 
the nature of randomness and is not only a function of 
the elastic properties. 



The anti-ferromagnetic Ising model on a rigid triangu- 
lar lattice is geometrically frustrated since the energy of 
the three bonds on each triangular plaquette of the lat- 
tice may not be simultaneously minimized [38l l39| . This 
leads to a highly degenerate ground-state and thus to un- 
conventional phases of matter [4fj, SH, S3, EH • Allowing 
the lattice to deform may relieve this frustration and lift 
the ground-state degeneracy 0, EE] ■ Recent experimen- 
tal j4(| and theoretical [47| studies have shown that frus- 
tration is only partially relieved and that such systems 
exhibit glassy behavior, dramatically slow down, and fall 
into metastable disordered configurations. In order to 
allow the system to obtain a disordered deformation, we 
follow [IH and assume a pair potential for Eq. ([9]) of the 
form: 

4>ij (rtj) = - J Oi<Jj [1 - e (n 3 - a)} + ]p {r l} - a) 2 , (43) 

Here we assume the spin variables tr, = ±1 are ran- 
dom. The magnetic interaction is taken to be anti- 
ferromagnetic J < 0, e > controls the magneto-elastic 
coupling strength, and 7 > is the stiffness of the uni- 
form springs connecting each nearest-neighbor pair. Note 
that some previous studies of the density of vibrational 
states of elastic networks (see for example [27]]) focused 
on harmonic lattices with random spring constants (a 
multidimensional version of the one-dimensional analysis 
provided in Section [rV]) . In the models we use here, the 
strength of the interaction is randomized by frustration 
and non-harmonic terms in the potential and therefore 
arises more naturally; we do not need any assumptions 
about the distributions that govern the interactions. In 
this sense the disorder in this model and its derivatives 
are similar in nature to the disorder in glass forming 
molecular systems. Another important difference with 
respect to conventional lattice models is the effect of off- 
lattice positional disorder on the terms that depend on 
the relative distance and orientation between two parti- 
cles in the matrix , which in the present case reads 
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FIG. 6: Color online: Typical realization of the elastic trian- 
gular anti-ferromagnet (|43p . Solid points represent up spins 
and empty points represent down spins, dashed lines connect 
interacting up spins and solid lines connect interacting down 
spins. Dotted lines connect pairs of interacting up and down 
spins. 



+ ■/gjgjc + 7(r <j --a) :r (M) 

We calculated the density of states for 20 realizations 
with 6400 particles each. Each realization of the system 
was initiated by positioning the particles on a triangular 
lattice with periodic boundary conditions. Each parti- 
cle was assigned a random spin value and the energy of 
the entire network was minimized, using the conjugate- 
gradient method (49[. The minimization was achieved 
by changing the coordinates of the particles, keeping 
the interaction between the original nearest-neighbors 
only, and keeping the spin values fixed. A typical re- 
sulting configuration is shown in Fig. [6] The density 
of states was then calculated using the eigenvalue Eq. 
([Ti?|) . The square-roots of the eigenvalues were collected 
in bins and the histogram recorded. To compute g{ui)/uj) 
the most precise method turned out to be calculating 
g(u)/u! = 2G(u! 2 ) where G{lo 2 ) is the histogram of the 
eigenvalues themselves. In order to compute g(u>)/tu at 
lu = we employed Eq. ([5]) and the elastic moduli com- 
puted below. The same method was used for all the 
models listed below. Throughout, we set 7 = 1, m = 1, 
J = 1, a = 1 and measure the density of states for var- 
ious values of e and of the other parameters defined for 
the subsequent models. 

Figure[7]shows the effect of disorder on the distribution 
of angles 9 between the inter-particle bonds and the x 
axis. This angle determine the value of the term Hij <8> 
riij. In a perfect triangular lattice these angles take six 
discrete values 9i = ir/3i, (0 < i < 5). In the disordered 
system these angles have a smooth distribution, and due 
to isotropy it is sufficient to consider the distribution of 
of the angles of one bond, say between [— tt/6, it/6]. 
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FIG. 7: Color online: Distribution of angles between nearest- 
neighbors in the elastic triangular anti-ferromagnet (|43|l for 
various values of e, see inset. 
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FIG. 8: Color online: Vibrational density of states in the 
elastic triangular anti-ferromagnet (|43[) for various values of 
e, see inset. 

Figure [H] shows the effect of disorder on the density of 
states. For the ordered triangular lattice the density of 
states exhibits van Hove singularities 50]. The most ob- 
vious effect of disorder is the smearing of the singularities 
and the flattening of the density of states. This results 
in filling the gaps between the singularities but also in 
some modes leaking to higher and lower frequencies. In 
particular, there is a change in the density of states at 
low-frequencies compared to the tail that characterizes 
the ordered lattice. It is important to note that when e 
becomes too large, the network begins to fold upon itself. 
In a more realistic model, say with next-nearest-neighbor 
interactions, where the particles are not physically linked 
to each other this folding is relieved by changing the co- 
ordination number (i.e. number of neighbors) . Below we 
will also study the effect of randomizing the coordination 
number. 

To emphasize the deviation from Debye's model we 
examine in Fig. [5] the density of states divided by the 
prediction of Debye's model, which for d — 2 is linear 
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FIG. 9: Color online: Vibrational density of states nor- 
malized by Debye's prediction in the elastic triangular anti- 
ferromagnet (|43[) . Symbols are the same as in Fig. [8] 

in frequency. A peak at low frequencies is observed and 
its position shifts to lower frequencies as the disorder in- 
creases. However, its height decreases upon increasing 
disorder. Thus in this model there is a negative correla- 
tion between the magnitude of disorder and the amount 
of the deviation from the Debye model (see also Figs. [2] 
and [3]) . Note that in this example it is hard to notice 
the deviation from the Debye model at low frequencies 
without dividing the density of states by to. 

Examining Figs. [5] and [5] we note that the density of 
states reaches zero at zero frequency in accordance with 
(j4|). Dividing by u> we observe a finite limit in Fig. [9] 
This behavior follows from Eq. J4J which predicts such 
a finite limit at d — 2. 




FIG. 10: Color online: Density of states for the model with 
non-linear elasticity (|45[) with k = 0.25 and different values 
of e, see inset. 




U! 

FIG. 11: Color online: Density of states for the model with 
non-linear elasticity (|45p (Fig. I10[) normalized by Debye's 
model. Symbols are the same as in Fig. [10] 



B. Non-Linear Springs 

Here we investigate the effect of random contributions 
to the harmonic part of the potential. This will bring us 
closer to generic systems. There is more than one way of 
doing so, and we therefore consider two different models 
for the inter-particle potential. The first has the form 

<f>ij(nj) = -J aw [1 - e (rij - a)} 

+ \l {Uj - a) 2 + \k (r tf - af (45) 
The harmonic term now reads: 

fl'tru) = 7 + Mnj - a) ■ (46) 

Due to the fluctuation in the inter-particle distances 
around a, this term fluctuates around an average value 
7- 

We repeated the procedure described above for calcu- 
lating the density of states for k = 0.25,0.5 and 1, and 
for various values of e. We observed the same qualitative 



behavior for all n values and present in Figs. [TOl andfTTI 
the raw density of states and the result after normaliz- 
ing by Debye's prediction for n = 0.25. As with the first 
model, we see excess modes at low frequencies 

C. Magneto-elastic coupling 

The second way to modify the elastic triangular anti- 
ferromagnet ()43|) is by adding a non-linear separation 
dependence to the magneto-elastic coupling term: 

<t>ij(nj) = -J o-.iO-j [1 - e (nj - a) + ~v{rij - a) 2 ] 

+ \l in 3 -a) 2 . (47) 
The harmonic term in this case reads: 

4>'l j {T ij )= 1 -Ja i a j y. (48) 
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FIG. 12: Color online: Effect of the non-linear magneto- 
elastic coupling (|47[) on the density of states for v — 0.3 and 
different e values, see inset. 



FIG. 14: Color online: Density of states of the diluted elastic 
triangular anti-ferromagnet (|49[) with e = 0.1 and various p 
values, as indicated in the inset. 




FIG. 13: Color online: The density of states with non-linear 
magneto-elastic coupling (Fig. I12[l divided by the Debye be- 
havior. The symbols are the same as in Fig. 1121 



The density of states for this case is shown in Figs. 
[T21 and IT51 for a representative value of v = 0.3. Quali- 
tatively similar results were obtained for v = 0.45. We 
recognize in these figures a smoothing of the Van-Hove 
singularities with redistribution toward both lower and 
higher frequencies. As before, we see a peak at low fre- 
quencies which moves toward lower frequencies when the 
disorder parameter e is increased. We will next examine 
the effect of topological disorder and see that this type 
of disorder has a much more pronounced effect on the 
density of the low-frequency states. 



the coordination number near isostati city for determining 
low- frequency vibrational modes [5l|, H2, HH . In order to 
account for these effects we use the elastic triangular anti- 
ferromagnet (T43|) but include the possibility of a missing 
link between two neighbors: 



,{-Jct^[1-£ (r y - a )] + - 7 {r l3 -a) 2 } (49) 



Where g%j is with probability p and 1 with probability 
1 — p. We use p close to to avoid rigidity percolation 
[541 155| and keep at least 3 bonds per particle in order to 
avoid floppy modes (modes of zero frequencies) . We thus 
create a sparse network with a local coordination number 
varying between 6 and 3. This model is a modification 
of the model described in [56[ which studied the effect 
of disconnecting links of a harmonic triangular lattice. 
In contrast to that model, our model introduces disorder 
in the equilibrium positions of the particles as well as 
in their coordination number. We solved this model as 
before, by first finding the lattice deformation that locally 
minimizes the mechanical energy. In this model the effect 
of disorder on the low-frequency domain of the density 
of states is much larger than before, as seen in Figs. [T4l 
and[l5l 

The slope of the density of states at low frequencies, 
although linear, as expected by Debye's theory, is very 
different from the slope of the density of states for the 
perfect lattice (see also Fig. @|. However this modified 
Debye point is consistent with the system's elastic mod- 
uli, as will be descried in the following section. 



D. Randomly diluted elastic triangular 
anti-ferromagnet 

It has recently been argued that the glass transition 
involves a change in the number of neighbors that each 
particle has [1, |g, 0, H, H, [Hi [H| • Moreover, recent work 
on jammed sphere packings has indicated the relevance of 



VI. ELASTIC MODULI 

To better understand the differences in types of ran- 
domness and their effect on the frequency redistribu- 
tion we consider here the elastic moduli of the two- 
dimensional models treated above. Contrary to the den- 
sity of states, the elastic moduli are global measures of 
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FIG. 15: Color online: Density of states of the diluted elastic FI G- 16: Color online: Shear moduli of the three models Eqs. 
triangular anti-ferromagnet g9| (Fig. normalized by the (SSI) (circles) , (45]) (squares) and (HTJ) (diamonds) for different 
Debye prediction. Symbols are the same as in Fig. M values of the disorder control parameter e. 



the response of the system to external mechanical pertur- 
bations. Nevertheless there exists an interesting relation 
between these global properties and the frequency redis- 
tribution. 

Measurements of the shear modulus were done by ap- 
plying affine shear transformations, minimizing the en- 
ergy after each step using the Lees-Edwards periodic 
boundary conditions, in order to measure the stress. Af- 
ter each minimization the stress for each particle was 
measured directly from its microscopic definition and the 
mean stress was computed as a sum over all particles. 
Next, the mean stress as a function of the strain was cal- 
culated, and the shear modulus was extracted from the 
numerical derivative. The bulk modulus was measured 
by decreasing the volume and measuring the diagonal 
part of the stress tensor (the pressure) . The elastic mod- 
uli were used to calculate the Debye point at lo = 0. 

We first measured the elastic moduli for the first three 
models Eqs. (g3]), (gSJ) and (g7|) for different values of e 
(see Figs. [TBI and IT7| . The results are somewhat unex- 
pected. In all three models the shear modulus increases 
when disorder is increased, while the bulk modulus de- 
creases. Thus in these three models we cannot say that 
the system softens or hardens, since one elastic modulus 
decreases while the other increases. 

For the diluted network both elastic moduli de- 
crease with increasing p, see Figs. [T8landfT9l This is very 
physical; cutting bonds must result in true softening of 
the system. Note that for a fixed value of p the qual- 
itative behavior with e is similar to the previous three 
models. 

We see that this last model differs from the previous 
three in having clear softening when the parameter p in- 
creases. One way to take into account both moduli in 
discussing the softening of the system is by focusing on 
the Debye point 



lim 



d-l 



(50) 
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FIG. 17: Color online: Bulk modulus of the three models 
Eqs. (J43j) (circles) , (|45} (squares) and (|47|) (diamonds) for 
different values of the disorder control parameter e. 
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FIG. 18: Shear modulus of the diluted model Eq. (|49)l as 
a function of p and e. The shear modulus increases when e 
increase but it decreases when p increases. 
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FIG. 19: Color online: Bulk modulus of the diluted model 
Eq. (|49[1 as a function of p and e, symbols as in Fig. 1181 The 
bulk modulus decreases when e increases but it decreases like 
the shear modulus when p increases. 
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FIG. 20: The Debye frequency of the three models Eqs. (l43|) 
(circles) , (|45|l (squares) and (I47p (diamonds) for different 
values of the disorder control parameter e. 

We computed the Debye frequency for the four models 
at hand, and the results are presented for the first three 
models in Fig. fSD] and for the fourth model in Fig. fJU 

We see from Fig. [501 that the Debye frequency is prac- 
tically constant for the first model, and slightly increases 
with e for the second and third models. For the fourth 
model (Fig. [3T]) the Debye frequency softens dramat- 
ically when the parameter p is changed. The disorder 
governed by the parameter e almost does not change the 
Debye frequency also in this model. 

VII. DISCUSSION AND CONCLUSIONS 

The main conclusion from the one-dimensional and 
two-dimensional examples treated above are as follows: 

1. Both ordered and amorphous solids exhibit peaks 
in their density of states. In ordered solids these 



2.9 




P 



FIG. 21: Color online: The Debye frequency of the diluted 
model Eq. (|49[) as a function of p and e, symbols as in Fig. 
EEH1 

peaks are understood as van Hove singularities. In 
amorphous solids these singularities are smoothed 
out, providing higher amplitudes to both lower and 
higher frequencies. 

2. For both the crystalline and the amorphous exam- 
ples the comparison with the Debye model shows 
complete agreement only at lu — > 0, independently 
of the dimensionality. Within the Debye model we 
expect g(uj)/ui d ~ l to be constant. This seems to be 
never the case. 

3. Dividing the computed density of states by Lo d ^ 1 
reveals the so-called "Boson peak". Its position 
and amplitude depend on many details; in one- 
dimensional cases we showed how it depends on the 
statistical distribution of the spring constants. In 
two-dimensions we showed how it depends mainly 
on the spatial disorder and on the coordination 
number, with the latter being dominant. In this 
sense there is nothing universal about the Boson 
peak. We cannot even tell a-priori whether increas- 
ing disorder might increase or decrease the ampli- 
tude of the Boson peak, cf. Fig. [22] 

4. We cannot discern any clear correlation between 
the Boson peak and the elastic moduli. In one di- 
mension we showed (Fig. [5]) that three models with 
identical bulk modulus exhibit completely different 
redistributions of frequencies. In two-dimensions 
we showed for the first three models that the bulk 
modulus decreased with disorder whereas the shear 
modulus increased, contrary to expectations. The 
change in the Debye frequency is small; never- 
theless we have completely different frequency re- 
distributions. In the fourth two-dimensional model 
we considered, both elastic moduli decrease simul- 
taneously and we indeed saw a pronounced redis- 
tribution to lower frequencies when the average co- 
ordination number changed. Again we see no sys- 
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FIG. 22: Color online: g(u})/u>) for the diluted model with 
different values of e. Note that the Boson peak amplitude 
is reduced when e which randomizes the particle positions is 
increased. 



tematic correlation with the behavior of the elastic 
moduli. 

Of course, all these conclusions concern the simple 
models discussed above. Nevertheless the phenomena 



discussed are not special to these models or even to amor- 
phous solids in general. An experimental connection be- 
tween shear modulus and the low- frequency behavior of 
the vibrational spectrum was given in [57| by analyz- 
ing the low-temperature specific heat. In contrast to 
the common view that excess in low-temperature specific 
heat (and, hence, in low-frequency modes) is special to 
disordered systems only, it was demonstrated that crys- 
tals and amorphous solids with almost the same shear 
modulus have quite different positions of the Boson peak. 
This is in accord with our conclusions that with the same 
Debye point we can have different redistributions of fre- 
quencies. 

In summary, it is quite possible that in a given family 
of amorphous materials, where the randomness is quite 
similar, there can be a correspondence between the re- 
distribution of frequencies and the shear modulus. How- 
ever this is not a general correlation, as we saw with 
the present examples. We saw that the actual density of 
states is a complicated function of many competing in- 
fluences. It is unlikely that one given parameter of what- 
ever nature (like the shear modulus) can capture this full 
complexity. The understanding of the density of states 
and its changes under modified interactions remains a 
theoretical calculation of significant difficulty. 
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